Superstatistical fluctuations in time series: Applications to share-price dynamics and 

turbulence 
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We report a general technique to study a given experimental time series with superstatistics. 
Crucial for the applicability of the superstatistics concept is the existence of a parameter j3 that 
fluctuates on a large time scale as compared to the other time scales of the complex system under 
consideration. The proposed method extracts the main superstatistical parameters out of a given 
data set and examines the validity of the superstatistical model assumptions. We test the method 
thoroughly with surrogate data sets. Then the applicability of the superstatistical approach is 
illustrated using real experimental data. We study two examples, velocity time series measured in 
turbulent Taylor-Couette flows and time series of log returns of the closing prices of some stock 
market indices. 

PACS numbers: 89.75.-k, 05.45.Tp, 89.65.Gh 



I. INTRODUCTION 

Superstatistical techniques as introduced in [l[ can be 
applied to large classes of nonequilibrium systems. The 
crucial assumption of this approach is that the dynam- 
ics of the system under study is a superposition of two 
dynamics with well separated time scales. The standard 
example is a Brownian particle moving through a slowly 
fluctuating environment (e.g. there is a slowly fluctuating 
temperature 1//?). The fast dynamics is then determined 
by the change of the velocity of the Brownian particle 
and the slow dynamics is determined by the change of 
the temperature of the environment. It is well known 
that the stationary velocity distribution of a Brownian 
particle in a constant environment is a Gaussian distri- 
bution with variance l/f3. The theory of superstatistics 
applies when the slow dynamics is so slow that the ve- 
locity distribution of the particle has time to relax to 
a Gaussian distribution between the changes of the en- 
vironment. As such, after a long time, the stationary 
velocity distribution of the particle is just a superposi- 
tion of Gaussian distributions weighted with a function 
/(/?). This /(/?) is the probability density to observe 
some value of (3. Depending on f((3), different results 
for the stationary velocity distribution will occur. An 
important question is which types of distributions will 
occur in 'typical' complex systems. In [2j, the authors 
give some probabilistic arguments in favor of three dis- 
tributions, the lognormal distribution, the gamma dis- 
tribution and the inverse gamma distribution. One can 
also derive these distributions from a maximum entropy 
principle [3]. The present paper will focus on techniques 
to extract the distribution /(/?) and the relevant super- 
statistical parameters out of real experimental data. 

Recent applications of superstatistical methods include 
hydrodynamic turbulence 0, 0, H, |(|, pattern form- 



ing systems 0], cosmic rays B, solar flares [9j], share- 
price fluctuations [l(| EH, Tl2l Il3ll. random matrix the- 
ory 0, random networks [16||. multiplicative- noise 
stochastic processes [l7[, wind velocity fluctuations [l8l |. 
hydro- climatic fluctuations [HI , the statistics of waiting 
times [13, and models of the metastatic cascade in 
cancerous systems 22]. For the present paper, we will 
concentrate on time series as generated by hydrodynamic 
turbulence and share-price fluctuations. 

The outline of the paper is as follows. In the next 
section we fix our notation. In section IIIII we discuss 
the techniques proposed in Q to extract superstatisti- 
cal parameters out of a given time series. We propose a 
modification to this method in section [TV] In section fVl 
we test our method using surrogate data sets and show 
that our modification improves on the previously pro- 
posed method [|J • In sections I VII and I VIII we apply the 
superstatistical approach to real experimental data from 
hydrodynamic turbulence and finance. An alternative 
procedure to study superstatistical time series is briefly 
examined in section IVIIII The final section contains a 
discussion of our results. 



II. NOTATION AND SUPERSTATISTICAL 
APPROACH 

The starting point is a given discrete time series u con- 
taining n data points. The different data points are de- 
noted as Ui with i = 1, 2, . . . n. The probability distribu- 
tion of the random variable u, extracted from the experi- 
mental data, is denoted as P(u). The total time series is 
divided into N equal slices of length A, with N = [n/ A\ 
where [x\ means rounding the value of x to the nearest 
lower integer. This implies that the Zth time slice contains 
the measurement points Ui with 1 + (I — 1)A < i < IA. 
We then define 'local' moments of order k as follows 
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with l = l,2,...N. (1) 
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In the following, we will assume that the first momentum 
of the total time scries u 
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u, 



(2) 



is equal to 0. If this is not the case we use the time series 
u — instead of u itself. 

Now we want to analyse the distribution P(u) using su- 
perstatistics. We assume that the time series u contains 
two different time scales r and T such that r/T << 1. 
Remember that we divided the time series u in different 
time slices of length A. The time scale r determines how 
fast 'local' equilibrium is reached in these different slices. 
The time scale T is the length of the slices. For every 
time slice, one can extract a 'local' probability distribu- 
tion for the variables with 1 + (I — 1)A < i < IA and 
A = T. When r/T << 1, it is a good approximation 
to assume that 'local' equilibrium is reached in a slice of 
length T. By this we mean that the 'local' distributions 
can be approximated by Gaussian distributions 



PT,l( 
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with @t,i — 1/(u 2 )t,;- Within this assumption, the dis- 
tribution P(u) is approximated by 



1 N 

P(u) ~p T {u) := ^^2ptA u ) 



(4) 



with N = \n/T\. In this way, one obtains N values for 
Pt,i ■ When N is large enough, one can replace expression 
Q by 



P(u) « p T (u) ps ptj{u) 




(5) 



with /t(/3) being the probability density that the value 
of the inverse variance in a randomly chosen time slice 
of length T equals f3. The obtained distribution Jt(P) 
depends on the long time scale T, because T determines 
the length of the time slices and as such the obtained 
values for the parameter (3. Therefore, the proper defi- 
nition of T is of crucial importance and a major part of 
the present paper will focus on this issue. 

In p|, a method was introduced to extract }t{0) out 
of a given experimental time series. The authors pro- 
posed to check the validity of their approach a posteriori 
by comparing the resulting distribution pxj{u), with the 
distribution P{u) extracted from the experimental data. 
Notice that the superstatistical approach includes two 
approximations. In the first step, one assumes the ex- 
istence of two time scales r and T such that in every 
time slice 'local' equilibrium is reached. In the second 
step one assumes the existence of a distribution /t(/3) 
replacing the summation in expression ^ by an integral. 



One would like to judge the validity of the first approx- 
imation before moving on to the second approximation. 
Therefore, we propose an extension of the method intro- 
duced in [2j to examine the quality of the first approx- 
imation. Afterwards, the validity of the second approx- 
imation can be checked by comparing the distributions 
Pt(u) and ptj(u). In the present paper we investigate 
the two steps of the superstatistical approach separately. 
We will test our arguments with surrogate data and with 
real experimental data sets. 



III. ORIGINAL DEFINITION OF TIME SCALES 

The correlation function of a time series u can be cal- 
culated as follows 
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The superstatistical short time scale r of the time series 
is defined by the exponential decay of C n ,t(u) @ 



C n , T (u) = e 1 C rh0 (u). 
In a function ka is introduced as 



KA = E K/ V' Witn KA .' 



(7) 
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Notice that ka,i is just the kurtosis of the Ith time slice. 
The superstatistical long time scale T is then defined by 
the condition 



(9) 



To understand the meaning of this definition, remember 
that the main assumption of superstatistics is the exis- 
tence of two well separated time scales. When this is 
true, one can approximate the distribution of the vari- 
ables in the Ith time slice by a Gaussian distribution (in 
this paper we always assume that 'local' equilibrium is 
associated with Gaussian behavior). When the variables 
are indeed locally Gaussian distributed with zero mean 
and variance 1//?a.z, the first four 'local' moments and 
ka.i are simply 

(u) A ,i = 0, (u 2 ) A ,i = (u 3 )a.i = 0, 

3 

(u 4 )a,1 = -fi2~, KA,! = 3. (10) 
<°A,Z 

The condition ([9]) implies that one is looking for a suit- 
able division of the total data set into time slices for 
which the variables are locally Gaussian distributed 0], 
with a variance that fluctuates from slice to slice. 

This definition will always result in a value for the 
long time scale whenever the kurtosis of the complete 
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time series K n> \ is larger than 3 @, [23[. This can be 
seen by considering two special cases. In the first case, 
A is so small that only one value of u is observed in 
each time slice. This results in K\ = Kij = 1 (with 
/ = 1,2, ...n). In the second case, A is so large that 
it includes the entire time series u. This results in n n = 
K n> i > 3. As a consequence, there exists a long time scale 
1 < T < n for which ka = 3 holds. With the definitions 
proposed in one could formally use superstatistics 
to analyse the distribution of any time series whenever 
t/T « 1 and K n> \ > 3. But clearly not every time 
series that fulfills these conditions contains time slices 
of which one can assume that the variables are locally 
Gaussian distributed. Therefore, in the next section we 
will derive an extra condition that should hold before we 
assume that a time series at hand can be described by 
superstatistics. 



IV. EXTRA CONSTRAINT 



@tj to the value of knt will be small as compared 
to the contribution of the term 3^-(u 2 )^j. There- 
fore, we expect the Gaussian approximation to hold when 
lei << 1 with 
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We showed that the value of the kurtosis of the complete 
data set n n ,i can be expressed as a function of the mo- 
ments of the times slices of length T only. The parameter 
e basically measures the contribution of the deviations 
from the Gaussian approximation in these time slices to 
the value of K n \. More details about the derivation and 
the interpretation of expression (|14[) can be found in ap- 
pendix [A] In the next section we will show that the extra 
constraint |e| « 1 resolves some of the ambiguities of 
the original approach [2|. 



Assume we found the value of A for which ka equals 3. 
This value is denoted as T. This means that the complete 
time series is divided in N = \n/T\ time slices of length 
T. The second and fourth momentum of these time slices 



(u 2 )tj 



and {u a )t,i respectively. Then one 

NT as 



(11) 



are denoted as 

can express the value of the function «a for A 
follows 




For large data series, the kurtosis of the complete data 
set K n< i is in good approximation equal to knt because 
N = [n/T\ w n/T. Formula HI]) shows that the value 
of K nt i can be calculated starting from the second and 
fourth momentum of the time slices of length T. Now we 
assume that these time slices contain variables that are 
approximately Gaussian distributed. We will determine 
the influence of the deviations from this approximation to 
the value of K n ,i. Remember that the moments of purely 
Gaussian distributed random variables are given by (|10[) . 
However, the variables of a time series are never perfectly 
Gaussian distributed. We define 9t.i as the deviation of 
the fourth momentum from 3(u 2 )^ l for the time slices of 
length T 

e T ,i = {u 4 )Tj-i{u 2 ) 2 T,i- (12) 

Then expression (fTTj) can be rewritten as follows 



knt 



« 2 >TJ + 0Tj] 



(13) 



When the Gaussian approximation is reasonable in the 
time slices of length T, the contribution of the term 



V. SURROGATE DATA 

The major problem of the definition © of the long 
time scale T is that it will always give a value for T when- 
ever the kurtosis of the complete time series K n> i is larger 
than 3. For example, a time series containing Gaussian 
distributed random variables of constant variance with 
just one outliner can have k„ i > 3 and t/T « 1 with- 
out being superstatistical. Also a time series containing 
q-Gaussian [24| distributed random variables with con- 
stant variance can be wrongly classified as being super- 
statistical. In this section we use surrogate time series to 
illustrate that the extra constraint |e| << 1 is a useful 
tool to decide whether a given time series is supersta- 
tistical. We study three examples, (A) data constructed 
by numerically integrating a discretized Langevin equa- 
tion, (B) time series containing outliners and (C) data 
sets containing q-Gaussian distributed random variables. 



A. Langevin-like surrogate data 

In order to construct realistic surrogate data sets, we 
numerically iterate a discretized Langevin equation 



U 



(15) 



where 7 is a constant, the Li correspond to Gaussian 
white noise with variance 1 and Pi are random variables 
with a certain distribution. This is the dynamics of the 
example discussed in the introduction, a Brownian par- 
ticle moving through a slowly fluctuating environment. 
When (3i = j3 is chosen to be a constant, one ends up with 
the standard Langevin equation. In this special case, it 
is known that the stationary probability density of u is 
a Gaussian distribution with inverse variance /3. The re- 
laxation time is I/7. We study the more general case in 
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FIG. 1: (Colour online) Extraction of the long time scale 
T from the condition ka = 3 for two surrogate data sets 
(see expression (|8]) for the definition of ka). The difference 
between the two time series is the width of the distribution 
of the fluctuating parameter f3, see Fig. [2] 



which f3i is a random variable that is gamma distributed 



fs(P) = 



r(a) 



(16) 



with a and 9 constants. However, in order to make the 
resulting surrogate time series superstatistical, the time 
scale T$ of the fluctuating parameter f3i must be larger 
than the time scale t$ — I/7 on which 'local' equilib- 
rium is reached. In this way, the probability distribu- 
tion of u will relax to a Gaussian distribution before the 
next change of the value of /Jj occurs. Therefore, we 
chose ts = 10 and Ts = 500. In practice we keep the 
value of Pi constant and iterate the Langevin equation 
over Ts = 500 steps. Then we change the value of (3i 
and iterate the discrete Langevin equation again over 
Ts — 500 steps. We repeat this procedure N$ = 500 
times. The result of this procedure is a time series with 
ns = TsNs — 250.000 data points. We construct two 
data sets, with different values of the parameters of the 
distribution $Tjo§, a = 5, 10 and 9 = 10/(a- 1). Then we 
apply the superstatistical approach to these time series. 

The resulting short time scales are r = 9.52, 9.21 for 
the values of the parameter a = 5, 10 respectively. Fig.Q] 
shows ka as a function of A for the two surrogate data 
sets. The black solid line and the red dashed-dotted line 
are obtained for a = 5, 10 respectively. The constraint 
ka = 3 results in the following values of the long time 
scale: T = 409, 521 for a = 5, 10. The corresponding 
values of |e| are equal to 0.014 and 0.004 respectively. 
Because both inequalities t/T « 1 and |e| << 1 hold, 
we conclude that the data series at hand can be de- 
scribed within the superstatistical approach, and that 
our method to extract the relevant time scales works very 
well. 
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FIG. 2: (Colour online) The black solid line and the red 
dashed-dotted line are representations of the distribution (|16(l 
with 8 = 10/(a — 1) and a = 5, 10 respectively. These dis- 
tributions are used to construct two surrogate time series. 
Applying the superstatistical approach to these data sets re- 
sults in an approximation /t(/3) for the imposed distribution 
/s(/3). The extracted distributions /t(/3) are represented by 
black dots and red squares for the data sets with a — 5, 10 
respectively. 



We continue by extracting /t(/3) out of the time series 
and compare the obtained distribution /t(/3) with the 
imposed distribution fs(/3)- The result of this calcula- 
tion can be seen in Fig. We conclude that there is an 
excellent agreement between fr(P) and fs(P) for both 
surrogate data sets. Finally, we construct the distribu- 
tion pt(u) and compare the obtained distribution pt{u) 
with the original distribution P(u), see Fig. [3J There is 
an excellent agreement between pt(u) and P(u) for both 
surrogate data sets. 
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FIG. 3: (Colour online) The black dots and the red squares 
represent the observed distributions P(u) for the surrogate 
data sets with a = 5, 10 respectively. The superstatistical 
approximations pr{u) are given by the black solid line and 
the red dashed-dotted line. The dashed lines are Gaussians 
with the same variance as P(u). 
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FIG. 4: Determination of the long time scale T from the 
condition ka = 3 for two surrogate data sets (see expression 
l|8jl for the definition of ka). The minimum of the solid line 
is reached for A = 150. The dashed line is a monotonic 
decreasing function of A. 



B. Outliners 

To illustrate the influence of an outliner to the super- 
statistical approach, we construct two surrogate data se- 
ries. The first time series contains ns = 25.000 Gaussian 
distributed random variables constructed by numerically 
iterating a discretized Langevin (|15[) with 7 = 1/5 and 
Pi = (3 — 1. Clearly, the extracted value of the kurtosis 
for this time series is approximately equal to 3. The sec- 
ond surrogate data series is constructed out of the first 
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FIG. 5: Plot of the obtained value of the superstatistical pa- 
rameter |e| as a function of the kurtosis of the used q-Gaussian 
surrogate data set. The solid line is a guide to the eye. 



by replacing one of the elements of the first time series 
by an outliner. As a consequence of this, the kurtosis of 
the second time series will generally be larger than 3 and 
is equal to 4.63 for our particular example. Therefore, 
expression ^ would formally give a value for the long 
time scale when applied to the second data series. This 
can be seen in Fig.[4]that shows a plot of ka as a function 
of A, computed from the first (dashed line) and the sec- 
ond (solid line) surrogate data set. For the first data set, 
l K A — 3|/3 is a monotonic decreasing function of A, while 
for the second data set the function \ka — 3|/3 reaches a 
single minimum at T = A = 150. We also computed the 
short time scale for this surrogate data set, r = 4.52. As 
a consequence, the inequality t/T « 1 holds. Finally 
we calculated the value of |e| for the second surrogate 
data series and obtained 0.41. Clearly, the inequality 
|e| << 1 does not hold in this case and the data series 
at hand cannot be described with the superstatistical ap- 
proach. After removing the outliner, one is left with a 
data series containing just Gaussian variables of constant 
variance (equivalent to the surrogate data set we started 
from) and no value for the long time scale can be found, 
see Fig. [4] This example shows that one has to be careful 
when using the superstatistical approach to study data 
series that contain outliners. Apparently, our criterion 
|e| << 1 can help to identify truly superstatistical dy- 
namics. 



C. q-Gaussian 

We continue by constructing surrogate data sets con- 
taining q-Gaussian random variables of constant vari- 
ance. This means that the random variables are dis- 
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tributed according to the following distribution 



P(u) 
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-1/(9-1) 



(17) 



with b and q constants. In the limit q — > 1, the dis- 
tribution (|17p approaches a Gaussian distribution with 
inverse variance b. The kurtosis of the distribution (|17p 
is 3(5 — 3q)/(7— 5q). Clearly, the kurtosis is larger than 3 
for q > 1. Therefore, expression ([9]) would formally give 
a value for the long time scale when applied to data sets 
containing q-Gaussian random variables with q > 1. Our 
aim here is to identify this time series as being not su- 
perstatistical, since there is no long time scale on which 
the variance changes. 

The kurtosis diverges in the limit q — » 7/5. That is the 
reason why the analysis presented in this section will be 
restricted to data series containing q-Gaussian random 
variables with 1 < q < 7/5. We use the generalized Box- 
Miiller method, proposed in [24[ , to construct the surro- 
gate data series. The number of elements of every data 
set is 715 = 25.000. The used values of q range from 1.05 
to 1.35, while the used values of b are 0.1, 1, 10. We calcu- 
late the value of A for which the constraint ka = 3 holds 
and determine the corresponding value of |e|. For ev- 
ery tuple (q, b), we repeated these calculations 500 times 
(we constructed 500 different data sets) and averaged the 
value of |e| over these different runs. The result of these 
calculations can be seen in Fig. [5l This figure shows a 
plot of the value of |e| as a function of the kurtosis of the 
used q-Gaussian distribution. The figure contains three 
curves, for three different values of b, 10 (•••), 1 (+ + +) 
and 0.1 (□□□). Clearly, the three curves are almost indis- 
tinguishable. More importantly, the curves are increasing 
functions of the kurtosis (determined by q). This shows 
that the extra constraint |e| << 1 will reject the time 
series containing purely q-Gaussian distributed random 
variables as being superstatistical when the value of the 
kurtosis is large enough. When the value of the kurto- 
sis is close to 3, the value of |e| will also become small. 
However notice that in these cases the corresponding q- 
Gaussian distribution is also better and better approx- 
imated by the Gaussian distribution. This means that 
the assumption of 'local' equilibrium can be relaxed to 
'global' equilibrium. As a consequence, there is no reason 
to the study these time series within the superstatistical 
approach, because the usual equilibrium statistical me- 
chanics (i.e. Gaussians with constant variance) can be 
used in very good approximation. 



VI. TURBULENCE DATA 

After testing our method with various types of sur- 
rogate data, we now apply our superstatistical analy- 




FIG. 6: Plot of the distribution P(u) (• • •) computed from 
a data series of velocity differences measured in turbulent 
Taylor-Couette flow. The superstatistical approximations 
Pt(u) and Ptj(u) are represented by the solid line and the 
dashed-dotted line respectively. 




FIG. 7: Plot of the extracted distribution /t(/3) (• • •) com- 
puted from a data series containing the velocity differences 
measured in turbulent Taylor-Couette flow. The solid line 
shows a fit of /t(/3) to a lognormal distribution, see expres- 
sion (I18J, with a = 0.5222 and 9 = -1.561. 



sis to real experimental data. We use time series ob- 
tained in an experiment performed by Lewis and Swin- 
ney (25j . They measured a single velocity component 
v(t) as a function of time t in turbulent Taylor-Couette 
flow for different Reynolds numbers Re. The stationary 
probability distribution P(u) of the velocity difference 
u(t) :— v(t + S) — v(t) at a given scale 8 is well-known to 
exhibit non-Gaussian behavior. 

The turbulence data sets we used from Swinney's ex- 
periment are not scale invariant. Contrary to other tur- 
bulence data there is no inertial range for these data and 
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no scaling exponents of velocity increments can be de- 
fined. But in it was shown that the existence of two 
effective time scales is well supported by the data. 

The values of the parameters for the data series used 
in the present paper are Re = 540.000, 5 = 16 and 
iik2x 10 7 . The two different time scales are extracted 
from the time series asr ~ 7.2 and T « 186, in agreement 
with the results of [2J. We also evaluated expression (TT4")) . 
The result is |e| = 0.0205 << 1. Hence we assume that 
the experimentally measured time series is superstatis- 
tical. We also constructed the distribution pr(u) which 
is an approximation for the real experimental distribu- 
tion P(u), see expression (|4]). Fig. [6] shows the excellent 
agreement between these two distributions. 

As mentioned in section II, the superstatistical mod- 
eling approach includes two approximations. Until now, 
we only discussed the first. After this approximation one 
ends up with N values of the 'local' inverse variance 0t,i- 
This series of numbers can be treated as a set of random 
variables with a certain distribution fr{P)- In the second 
step of the superstatistical modeling process one tries to 
find the best fit to fr(P) with some well-known distribu- 
tions such as, e.g. the lognormal distribution, the gamma 
distribution or the inverse gamma distribution, thus pro- 
ceeding to an analytic model. Fig. [7] shows our extracted 
distribution Jt(P)- A good fit to the data is obtained for 
a lognormal distribution 



M/3) 
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(18) 



with a = 0.5222 and 6 = -1.561 (see also Q). Finally 
we construct ptj(u) (see expression ([5])) and compare 
this distribution with pt(u). Fig. [6] includes plots of 
both ptj{u) and pt(u). To summarize, Fig. [7] shows 
the good fit of a lognormal distribution to /t(/3), while 
Fig.[B]shows the excellent agreement between pxj (u) and 
Pt{u). Both figures validate that lognormal superstatis- 
tics is indeed a good modeling approach to the data. 
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1 
2 
4 


<1 17 1.11 35.6 
1.28 33 1.12 22.6 
2.46 60 1.02 15.4 


<1 18 1.50 49.3 
1.27 32 1.28 28.0 
2.47 64 1.36 18.6 


1 
2 
4 


<1 17 0.02 6.6 
1.29 32 0.01 6.1 
2.45 73 0.05 5.6 


<1 19 0.01 6.3 
1.29 34 0.01 6.0 
2.45 70 0.10 5.5 



TABLE I: The values of the superstatistical parameters ex- 
tracted from data series containing the normalized log returns 
Ui of the DJI index and the S&P 500 index for three differ- 
ent values of 8. The used data sets are (top) unaltered m, 
(bottom) largest values of \ui \ removed. 
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FIG. 8: Plot of the normalized log returns m as a function 
of time i for the DJI index with 5 = 1. The plot shows just 
a small section (150 data points) of the complete data series 
(approximately 15.000 data points). Notice the outliner near 
i = 9500 which corresponds to the crash of the stock markets 
on October 19, 1987. 



VII. ECONOMICAL DATA SERIES 

In this section we study economical time series. The 
data sets are the daily closing prices Xi of the Dow Jones 
Industrial Average index (DJI) and the S&P 500 index 
for the period March 1950 to September 2008. This 
means that the total number of data points is of the 
order n ~ 15.000. In the literature one usually studies 
the statistics of the log returns Xi :— hi(xi+s/xi) with 
S = 1,2,... of the closing prices instead of the closing 
prices itself. We consider the normalized log returns 

Ui := (Xi (X)) (7<A^) - (X) 2 ) , (19) 

which have been rescaled to have variance 1. An example 
of a part of such data series with 5 — 1 for the DJI 
index can be seen in Fig. [5] This figure illustrates that 



sometimes outliners (crashes of the stock markets) are 
present in the data. 

In the first step of the superstatistical analysis, we cal- 
culate r, T and |e|. Additionally, we also calculate for 
every set the total kurtosis Ki >n . We performed these 
calculations for the DJI index and the S&P 500 index for 
three different values of 6 = 1, 2, 4, see Table Q] For every 
value of 5, the condition r/T « 1 holds, while the con- 
dition |e| << 1 does not hold. Notice also that the values 
of the total kurtosis strongly deviate from 3. We illus- 
trated in section fV Bl that outliners can strongly influence 
the results of the superstatistical analysis. Therefore, we 
continue by studying data sets in which the data points 
with the largest values of \m\ are removed. Indeed, it 
is well-known in mathematical finance that large jump 
events of prices play a special role and lead to modifi- 
cations of the Black-Scholes theory 26]. These events 
are often created by special circumstances of the market 
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FIG. 9: (Colour online) Plot of the distributions P(u) for the 
DJI index and the S&P 500 index for two different values of 8 
[S&P 500 index with 5 = 1 (•••) and S = 2 (□□□), DJI index 
with 8 = 1 (> > >) and 8 = 2 (ooo)]. The superstatistical 
approximations pr(u) are represented by the solid lines. The 
figure also shows the distributions pr,f for 3 different choices 
of /t(/9), the gamma distribution (dotted lines), the lognor- 
mal distribution (dashed dotted lines) and the inverse gamma 
distribution (dashed lines). The values of the corresponding 
parameters a and 8 are listed in Table El 



and are not explainable by superstatistics. We repeated 
the calculation of r, T and |e| for data sets excluding 
the outliners. In practice we removed data points with 
|ttj| > 7 (approximately 0.05% of the data points). The 
bottom part of Table U shows the newly obtained values 
for r, T, |e| and the total kurtosis K,\ >n when the outliners 
are removed. By comparing the values of the supersta- 
tistical parameters in Table [J one observes that r and 
T are hardly influenced by the outliners. However, no- 
tice the large decrease of the values of ki.„ and |e|, in 
particular for the data sets with S = 1,2. Removing ap- 
proximately 0.05% of the data points results in a decrease 
of the values of «i jn and |e| by a factor of the order ~ 5 
and ~ 100 respectively. We conclude that both condi- 
tions t/T << 1 and |e| << 1 hold for the economical 
data sets with S = 1,2 once the outliners are removed. 
Therefore, we will focus our analysis in the following part 
to these data sets. 

After the calculation of T one can construct the distri- 
butions pt(u) and /t(/3). Fig.[5]shows the original distri- 
butions P(u) of the normalized log returns together with 
the result pr(u) of the superstatistical approach after the 
first approximation. There is an excellent agreement be- 
tween P(u) and Pt(u) for all four cases. One also ends up 
with N w 15.000/T different values for the inverse vari- 
ance Pt,i- We constructed a histogram with these values 
and tried to approximate this histogram with some well- 
known distributions, the gamma distribution (j!6j) . the 
lognormal distribution (|18p and the inverse gamma dis- 
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FIG. 10: (Colour online) Plot of the extracted distributions 
/t(/3) for the DJI index and the S&P 500 index for two dif- 
ferent values of 8 = 1, 2 [S&P 500 index with 8 = 1 (•••) and 
8 = 2 (□□□), DJI index with 5 = 1 (>») and 8 = 2 (ooo)]. For 
all 4 distributions an approximation to the data in terms of 
a gamma distribution (dotted lines) , a lognormal distribution 
(dashed dotted lines) a inverse gamma distribution (dashed 
lines) is also shown. The values of the corresponding param- 
eters a and 8 are listed in Table HTl 



tribution 



MP) 



r(a) 



-6/0 



(20) 



These three distributions were motivated in Q. The re- 
sulting distributions are shown in Fig. [TUl We also con- 
structed the distributions ptj, see Fig.[S] By inspecting 
Fig. [9] and Fig. [10] we conclude that certainly the inverse 
gamma distribution is not a good candidate to represent 
/t(/3). However, it is much harder to distinguish between 
the lognormal distribution and the gamma distribution. 
Both distributions are a reasonably good approximation 
for /t(P) and result in a good agreement between pj-jiu) 
and pt{u). We conclude that the data sets under study 
can be described equally well by gamma superstatistics 
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S&P 500 


DJI 






G L I 


G L I 


1 
2 


a 

e 

a 



2.02 0.78 1.74 
1.05 0.56 2.45 
2.04 0.80 1.54 
0.98 0.49 2.01 


2.05 0.84 1.28 
0.96 0.51 1.64 
2.10 0.77 1.47 
0.88 0.44 1.83 



TABLE II: The values of the parameters a and 8, used to con- 
struct the distributions shown in Fig. [9] Fig. [10] and Fig. [11] 
The letters G, L and I stand for the gamma distribution (|16|l . 
the lognormal distribution (|18[1 and the inverse gamma dis- 
tribution (|20|) . respectively. 





S&P 500 


DJI 


<5 


r T \e\ ki,„ 


T T £ Kl, n 


1 
2 
4 


<1 8 0.30 6.6 
<1 9 0.30 6.1 
<1 10 0.30 5.6 


<1 9 0.32 6.3 
<1 10 0.32 6.0 
<1 11 0.30 5.5 



TABLE III: The values of the superstatistical parameters ex- 
tracted from data series containing the randomly shuffled nor- 
malized log returns m of the DJI index and the S&P 500 index 
for three different values of S (with largest values of \m\ re- 
moved) . 



and lognormal superstatistics. 

We also studied the effect of a random shuffling of the 
data sets. A shuffling of the data keeps the kurtosis and 
the distribution P(u) unchanged but destroys the corre- 
lations in the time series. As a consequence, we expect a 
considerable change in the values of r, T and |e|. We ap- 
plied the shuffling operation to the data of the DJI index 
and the S&P 500 index for S = 1,2, 4, again with outlin- 
ers removed. We shuffled every data set 500 times and 
averaged the values of T and |e| over these different runs. 
The result of these calculations can be seen in Table IIIII 
As expected, the random shuffling operation causes a de- 
crease of the value of T and an increase of the value of 
|e| (compare Table HI and Table IIII[) . These observations 
confirm the fact that the original economical time series 
is superstatistical whereas the shuffled one is not. Let us 
connect this with the results obtained in scction lVCl We 
previously provided evidence that the distribution P(u) 
of the economical data sets can be approximated very 
well by ptj(u), see expression ([S]) with fr{P) being a 
gamma distribution. In this case, the integral of expres- 
sion ([5]) can be evaluated analytically as 

This distribution is a q-Gaussian distribution, see expres- 
sion fP?|). with 

a = — and 6 = b(q-l). (22) 

q — 1 2 

In section IV CI we studied the dependence of the value 
of |e| as a function of the kurtosis for surrogate data 
sets containing q-Gaussian random variables. Because 
the economical data sets are in good approximation q- 
Gaussian distributed, it is no surprise that the values of 
the kurtosis (~ 6) and |e| (~ 0.3) of Table [m] are of 
the same order of magnitude as the values that can be 
read off the curve shown in Fig. [5] In other words, the 
randomly shuffled log returns can be modeled in good ap- 
proximation with the q-Gaussian random variables stud- 
ied in section IVCl 

Finally, we studied how our superstatistical techniques 
were influenced by extreme events (outliners) such as 
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FIG. 11: (Colour online) Same picture as Fig. [9] but with 
outliners included. 



stock market crashes. We formally repeated our eval- 
uation of the distribution pt{u) for the data series in- 
cluding the outliners, see Fig. QTJ This figure still shows 
excellent agreement between the superstatistical approxi- 
mation pt(u) and the original distribution P(u) for small 
and intermediate values of u but fails to reproduce the 
data in the region \u\ > 10. The reason for that is quite 
clear: The basic assumption of superstatistics, namely 
'local' equilibrium, is violated during periods of very high 
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FIG. 12: (Colour online) Difference between the two defini- 
tions of the long time scale T for time series containing the 
normalized log returns (with outliners removed) of the S&P 
500 index for 4 different values of 5. The average kurtosis ka 
of these time series is plotted for 5=1 (>>>), 5 = 2 (□□□), 
5 = 3 (•••), 5 = 4 (ooo). One definition fixes the value of T 
by kt = 3 (dashed line), while the other definition uses the 
equation k t = 3T/(2 + T) (solid line). 



market volatility. Again we approximated the extracted 
/t(/3) with the 3 aforementioned distributions and con- 
struct the corresponding distributions prj{u). When 
outliners are included, as Fig. [TT1 shows, the best agree- 
ment between pt{u) and pxj is obtained for the gamma 
distribution /t(/?)- This analysis suggests that without 
eliminating the outliners gamma superstatistics is best 
suited to model realistic share-price dynamics. 

VIII. ALTERNATIVE APPROACHES 

The definition of a long superstatistical time scale T is 
not unique and alternative approaches can be considered 
in principle. In this section we introduce an alternative 
definition for the long time scale and briefly repeat the 
study of the time series examined in sections IVT1 and IVlTl 
thus testing for the robustness of our methods. 

Since superstatistics is based on a superposition of 
Gaussian distributions, one might consider a suitable 
weighting of relevant observables with local Gaussian 
statistics. An interesting alternative approach is the com- 
parison of the average kurtosis ka of the time series at 
hand with the expectation of the kurtosis k(A) of a time 
series with length A averaged over Gaussian statistics. In 
appendix IBl we show that k(A) = 3 A/ (2 + A). The long 
time scale T can then be extracted from the condition 

k t = «(T) = 3^;. (23) 

This definition for T coincides with expression ^ in the 
limit T — > oo. As a consequence, it is no surprise that for 



the turbulence data the resulting values of the long time 
scales almost coincide: T — 186 (kt — 3) and T = 170 
(kt — k(T)). We also checked that there is no visible 
change of the distributions /t(/3) ° r Pt(u) in Fig. [5] and 
Fig. [7] with this new method. 

For the economical time series the typical values of 
T are smaller as compared to the turbulence data, and 
hence stronger differences might arise in this case. Fig.[T2l 
shows the average kurtosis ka for time series containing 
the normalized log returns (with outliners removed) of 
the S&P 500 index for different values of 5 = 1,2,3,4. 
As expected, the obtained values of the long time scales 
are different for the two methods, T = /; 17; 36; 52 
(k t = k(T)) and T = 17; 32; 52; 73 (k t = 3) for the time 
series with S = 1,2,3,4 respectively. Here the slash / in- 
dicates that no solution is found, due to a lack of intersec- 
tion of the relevant curves in Fig. 12. Most importantly, 
however, we checked that again the corresponding distri- 
butions /t(/3) or pt{u) obtained for the new method are 
almost indistinguishable from Fig. [5] and Fig. [TU1 This 
shows that our results are robust and that our main con- 
clusions about the time series at hand do not change. 

In section ITVl we derived a constraint |e| << 1 that 
should hold if a given time series is superstatistical. For 
the alternative rule, expressions (fT2"l) and (fTl)) must be 
replaced by 

Ot.i = (u 4 )TJ-S-^j;(u 2 )h 

e - ii±I (24) 

3 T zU u2 )h' 

With this expression for e, we obtain |e| = 
/; 0.03; 0.003; 0.06 (k t = k(T)) for the time scries 
with S = 1,2,3,4 respectively, whereas previously we 
had with expressions (fT2|) and (fT4|) the values |e| — 
0.02; 0.01; 0.03; 0.05 (k t = 3). All these numbers are 
small, so our conclusions don't change and are robust, 
though it is noticed that the new method apparently fails 
to provide an e and T- value for the case 5=1. Further 
alternative definitions can be given, and the methods can 
be extended to include tests of statistical significance. 

IX. DISCUSSION 

In this paper we have presented a general technique to 
study a given experimental time series with superstatis- 
tics. The crucial assumption of superstatistical model- 
ing is that the dynamical description can be split into 
two levels that have a large time scale separation. Then, 
the complete data set can be divided into different time 
slices of length T for which the variables are Gaussian dis- 
tributed with a certain fixed inverse variance (3. However, 
the value of /? varies from slice to slice. If the number of 
time slices is large enough one can construct a histogram 
with these values and approximate this histogram with a 
known distribution. We proposed a general procedure to 
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introduce the different approximations of the supersta- 
tistical approach. Our method is as follows 

• Extract the short r time scale ([7]) , the long T time 
scale (|8I9|) and the value of e ([T4|) out of the time 
series at hand. 

• Check whether the following two inequalities 
t/T « 1 and |e| << 1 are satisfied. Then con- 
struct the distribution pt{u) given by (J4j) which is 
a first approximation for the original distribution 
P(u). 

• Search for a good fit to the histogram of the val- 
ues of (3 with a known distribution Then 
construct the distribution ptj{u) ([5]) which is a 
second approximation for P(u). Only when a good 
fit is obtained together with pxj(u) w pr{u) one 
can conclude that the second approximation of the 
superstatistical modeling approach is valid. 

We tested this method with several surrogate data sets 
and showed that our method is able to extract the cor- 
rect information out of a given data set. We then ap- 
plied the proposed techniques to two real experimental 
data series, velocity time series measured in turbulent 
Taylor- Couette flow and time series containing the nor- 
malized log returns of the closing prices of some stock 
market indices. For the turbulence data the inequalities 
t/T « 1 and |e| << 1 were immediately seen to hold, 
whereas for the share prices the outliners had to be re- 
moved first in order to obtain a realistic superstatistical 
model. We conclude that the superstatistical approach 
can be successfully used to study both data sets. 

Since the early work of Black and Scholes [27], vari- 
ous techniques borrowed from the field of the theoretical 
physics were successfully used to study the evolution of 
stock markets prices and their derivatives. Some recent 
work in the context of option pricing involves for example 
the use of continuous-time random walks [HI, [2§| , pertur- 
bation expansions around the Black-Scholes formula (30j 
and the use of path integrals [3l|. Other authors focus 
on stochastic volat ility and its extraction from a long se- 
quence of data [10l.l32l.l33l]. Common in all these papers is 
the observation that the volatility of the log returns Ui of 
the stock market prices is a stochastic variable with cer- 
tain distribution. However the type of this distribution 
is still under debate. Often the volatility is defined as 
the average of \ui\ over a time slice with certain length 
A. This results in N = n/A values for the volatility, 
with n the total number of data points. Then one con- 
structs a histogram with these N values and searches for 
the best fit to this distribution using some known distri- 
bution. Usually, one examines different values of A and 
observes that the result of the fitting procedure is not 
crucially dependent on the arbitrary choice of A. Notice 
the differences with the method presented in this paper. 
The length of our time slices coincides with long time 
scale T. Therefore, in the context of superstatistics, the 
value of A is fixed by the definition of T. We also use 



another measure for the volatility. Instead of studying 
the average of |itj| over a time slice we study the inverse 
of the average of uf . The reason for this is that in the 
second step of our method we need the distribution of the 
latter variable to construct the distribution ptj{u). Our 
careful analysis of the different approximations of the su- 
perstatistical approach shows that for the hydrodynamic 
turbulence data our techniques can be applied directly 
whereas for the economical time series it is better to first 
remove the outliners. When outliners are removed, our 
data sets can be described equally well by gamma super- 
statistics and lognormal superstatistics. When outliners 
are included, gamma superstatistics seems to do the best 
job. 

Our approach is inspired by the theory of statistical 
hypothesis testing. In this context, the null-hypothesis 
is the assumption that the data series at hand can be 
described by the superstatistical approach. Then one 
has to calculate the values for r, T and e. One accepts 
the null- hypothesis when T exists and t/T « 1 and 
|e| << 1. An interesting topic for future research is a 
statistical analysis of the threshold behavior of t/T and 
|e|. Also other conjectures to falsify the null-hypothesis 
can be examined. Generally this work will further help 
to understand the behavior of complex systems with time 
scale separation, making direct contact with experimen- 
tal measurements. 



Acknowledgments 

We are very grateful to Harry Swinney for providing us 
with the measured turbulence data. E.V.d.S. gratefully 
acknowledges private discussions with Bart Partoens. 



APPENDIX A: 

The superstatistical long time scale T is defined as the 
value of A for which ka equals 3, see expression for 
the definition of ka- Assume that we have found the 
value of the long time scale. This means that the com- 
plete time series is divided in N — \n/T\ time slices of 
length T. The superstatistical model assumption is valid 
when the variables within these time slices are in good 
approximation Gaussian distributed with a certain vari- 
ance (u 2 )t,i, with I — 1...N. We defined 9t,i as the 
deviation of the fourth momentum from 3(w 2 )|i ; in the 
lth time slice of length T, see expression (IT21 . Continue 
by dividing the complete time series in \ n/sT\ time slices 
of length A = sT with s = 2, ... N. We will show in this 
appendix that the value of k s t can be expressed as a 
function of (u 2 }t,i and 0t,i only. The contribution of Ot.i 
to the value of k s t is a measure for the deviations from 
the Gaussian distribution in the time slices of length T. 

First notice that the moments within the lth time slice 
of length sT can be written as a sum of the moments 
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within the time slices of length T 

isT 



{u k )sT,i 



1 g 

i=l+(Z-l)sT 
1 5 

= ~ D< ufc >T,.(I-l) 



(Al) 



i=i 



with ^ = 1,2,... [n/sT\ . This gives for the fourth mo- 
mentum and K S T i 



1 s 

(u 4 )sT.l = -^<U 4 ) T)S (/-l)+j 



^sT,Z - 3 



*(/-l)+3 ^2 



T,a(J-l)+i. 



2\2 {BsT,l + ©sT.i) 



(u 2 > 



(A2) 



sT.Z 



where 



sTJ 



S l^X ^T,s{l-X)+j \ S 7~i Pt,s{1-X)+3 



1 S 



(A3) 



To understand the meaning of -B s t.z and 0st,z, we study 
an example. Assume I — 1 which means that we look for 
the first time slice of length sT. This time slice contains 
the first s time slices of length T, with corresponding 
values of the inverse variance /3t,i, ■ ■ ■ , Pt,s respectively. 
For this special case, the expressions for B s t : i and Q s t,i 
simplify to 



Fig. Q] For both curves, the contribution of J2i ©sT,z 
is small compared to the contribution of B s t.i and 
can be ignored (remember these graphs are obtained for 
surrogate data, where we imposed Gaussian distributed 
random variables). However, the values of -B s t,( are 
clearly different for the two graphs, because of the differ- 
ence in the imposed distribution fs{f3), see Fig. [2j The 
narrower the distribution fs(P), the smaller the fluctua- 
tions of the parameter ftrj , the less profound the mini- 
mum of the graph of \ka — 3 j /3. 

The contribution of & s t,i to (k 8 ti — 3)/3 vanishes 
when the Gaussian approximation in the time slices of 
length T holds. So this term measures the deviations 
from pure Gaussian behavior in these time slices. We 
propose to study the difference between the exact value 
of (k s t — 3)/3 with s = 2, 3, . . . and expression (|A5[) with 
©sT,; = for all I to decide whether a given time series 
can be described within the superstatistical approach. 
Clearly, this difference will show fluctuations as a func- 
tion of s. Therefore, it is reasonable to evaluate this dif- 
ference in the limit s — ► N = \ n/T\ in order to quantify 
the influence of @ s t,i to the exact value of {k s t — 3)/3. 
For the special case s = N, formula (jAljl simplifies to 



{u k )nt,x = 



(A6) 



Then we obtain for knt.i (or knt, because there is only 
one time slice for this special case) 



K NT 



[Z{u 2 ) 2 T,j + Btj] 



(A7) 




(A4) 



B S T,z is the variance of 1 /pr,j calculated over s values of 
this parameter. s t,z is the average of 0t,j over s values 
of this parameter. B s t.i vanishes when the fluctuations 
of @T,j are small. O s t,z vanishes when the Gaussian ap- 
proximation for the time slices of length T holds. 

We continue by summing the expression (k s t.i — 3)/3 
over all the time slices \_n/sT\ . This results in an expres- 
sion for (k s t — 3)/3 

[n/sT\ 



~~ ^ = 1 TW\ E 7~2\2 ( B *T,l + @sT,i) ■ 

[n/sT\ ^ (w 2 )st,z 

(A5) 



Notice that this formula explains the difference in pro- 
foundness of the minimum between the two graphs of 



Notice that for large data series, the kurtosis of the com- 
plete data set K n x is in good approximation equal to 
k nt because N = \n/T\ n/T. When the Gaussian 
approximation is reasonable in the time slices of length 
T, the contribution of the term ^ . 8tj to the value of 
knt will be small as compared to the contribution of the 
term 3J2j( u2 )Tj- Therefore, we expect the Gaussian 
approximation to hold when |e| << 1 with 
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"r,j 



3 EL 



(A8) 



The value of e basically measures the contribution of the 
deviations from the Gaussian approximation in the time 
slices of length T to the value of the kurtosis of the com- 
plete data set. 

Finally we illustrate the physical interpretation of the 
formulas we derived above by applying them to the real 
experimental data that we already studied in section [VTl 
We repeated the calculation of \k& — 3|/3 as a function 
of A with formula © . The result can be seen in Fig. [T3] 
(solid line). Then, we evaluated expression (| A5[) with 
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FIG. 13: The solid line shows the value of \ka — 3]/3 calcu- 
lated with expression (|8}. The evaluation of expression (|A5|) 
with O a T,i = for all I is represented by (+ + +). The dif- 
ference between the solid line and the curve represented by 
(+ + +) is a consequence of the deviations from the Gaussian 
approximation in the time slices of length T. 



A averaged over a Gaussian statistics 



k(A) 



2tt 



+ 00 



+ CO 



■ exp 




dxA x 



(Bl) 



To calculate these A integrals it is convenient to consider 
a transformation to hyperspherical coordinates [3^ | 



X\ = r cos (pi 

x-2 = r sin <f>\ cos <f>2 

xa-1 — rsin^i . . .sin^A-2 cos^a-i 



(B2) 



xa = r sin 0i . . . sin </>a-2 sm 
with r £ [O..oo], <pi, ■ ■ ■ ,(f>A-2 S [0..ir] and 0A-1 <= 
[0..27r]. The Jacobian of this transformation is 



©sT,i = for all I (+ + +). Clearly, the difference be- 
tween the solid line and the curve represented by ( + + + ) 
is small. This shows that the contribution of the term 
J2i ®T,j / (u 2 )t j t° the value of k s t is small as compared 
to the contribution of the term J2j Bt.j/ ( m2 )tj - Hence 
the Gaussian approximation in the time slices of length 



T is reasonable because the term J2i ®T,j/(u")T j mea- 
sures the deviation from this approximation. The value 
of |e| corresponds to the difference between the solid line 
and the curve represented by (+ + +) in the limit s — > N. 
For this example we obtain |e| = 0.0205 << 1. 



r A - x sin(0 1 ) A -^ sin(0 2 ) A - J . . . sin(0 A _ 2 ). (B3) 

With the use of the well-known formula for the volume 
of a A-dimensional sphere with radius R 



r A/2 



r(A/2 + 1) 



(B4) 



one can easily check that the expression for k(A) simpli- 
fies to 



APPENDIX B: 

In this appendix we calculate the mathematical expec- 
tation of the kurtosis k(A) of a time series with length 



k(A) = 3 



A 



2 + A 



(B5) 
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